biosignalsnotebooks | download button biosignalsnotebooks | binder server biosignalsnotebooks | home button biosignalsnotebooks | contacts button biosignalsnotebooks | github button

Digital Filtering - A Fundamental Pre-Processing Step
Difficulty Level:
Tags pre-process☁filter

The acquired electrophysiological signals have always two intrinsic components. The signal we really want to acquire/study and noise, i.e. the acquisition component that is not relevant for the study purposes.

Noise can have different origins, such as in random events or due to voluntary/involuntary movements of the subject under analysis that affect the signal acquisition .

So, filtering is a fundamental step that needs to be applied to the signal, in order to ensure the maximisation of signal to noise ratio . Filtering can be achieved by hardware, having the analogical systems a great relevance, or by software using digital filters.

In this Jupyter Notebook it will be demonstrated how to digital filter the signal.

1 - Importation of the needed packages

In [1]: 
# biosignalsnotebooks own package for loading and plotting the acquired data
import biosignalsnotebooks as bsnb

# Scientific package
from numpy import array, mean, average, linspace, where
from numpy.random import normal
In [2]: 
# Hide warnings.
# https://groups.google.com/a/continuum.io/forum/#!topic/bokeh/h817iNS2twk
from IPython.display import HTML
HTML('''<script>
code_show_err=false; 
function code_toggle_err() {
 if (code_show_err){
 $('div.output_stderr').hide();
 } else {
 $('div.output_stderr').show();
 }
 code_show_err = !code_show_err
} 
$( document ).ready(code_toggle_err);
</script>
To toggle on/off output_stderr, click <a href="javascript:code_toggle_err()">here</a>.''')
Out[2]:
To toggle on/off output_stderr, click here .

2 - Load of acquired ECG data

In [3]: 
# Load of data
data, header = bsnb.load_signal("ecg_4000_Hz", get_header=True)

In the following cell, some relevant information is stored inside variables. This relevant information includes the mac-address of the device, channel number and signal acquisition parameters such as resolution and sampling rate.

For a detailed explanation of how to access this info, the "Signal Loading - Working with File Header" Notebook should be consulted.

In [4]: 
ch = "CH1" # Channel
sr = header["sampling rate"] # Sampling rate
resolution = header["resolution"] # Resolution (number of available bits)

3 - Generation of signal power spectrum by Fast Fourier Transform (FFT)

With this step is possible to observe the frequency composition of ECG signal.

3.1 - Store the desired physiological data (channel 1) int an individual variable

In [5]: 
# Acquired data
signal = data[ch]

3.2 - Removal of continuous component from our signal (baseline shift through the subtraction of the average value)

This task ensures more stability of our filtering system.

In [6]: 
# Baseline shift.
signal = array(signal) - mean(signal)

3.3 - Generation of the power spectrum

In [7]: 
# Power spectrum
freq_axis_1, power_spect_1 = bsnb.plotfft(signal, sr)

4 - The informational content of ECG signal is typically contained below the 40 Hz frequency component

With the next representation we can conclude that exist some unwanted information out of this frequency band.

In [8]: 
fig_1 = bsnb.plot_informational_band(freq_axis_1, power_spect_1, signal, sr, band_begin=0.5, band_end=40, legend="ECG Power Spectrum", 
                                     x_lim=[0, 100], y_lim=[0, 5e6], show_plot=True)

5 - Application of a low-pass filter in order to be excluded the unwanted information above the 40 Hz frequency component

Some low-frequency noise can be present at [0, 0.5] Hz frequency band. To exclude it we can follow an identical procedure, but, instead of applying a low-pass filter, it should be used a band-pass filter for the frequencies inside [0.5, 40] Hz.
For now, we focused on the more problematic type of noise, i.e., the high frequency noise

In [9]: 
# Digital lowpass filtering with a cutoff frequency f of 40 Hz
filter_signal_1 = bsnb.lowpass(signal, f=40, order=1, fs=sr)

# Power spectrum
freq_axis_2, power_spect_2 = bsnb.plotfft(filter_signal_1, sr)

6 - Comparison of the power spectrum of original and filtered signal

In [10]: 
bsnb.plot_before_after_filter(signal, sr, band_begin=0.5, band_end=40, x_lim=[0, 100], y_lim=[0, 5e6], show_plot=True)
Out[10]: 
[[Figure(id='1180', ...), Figure(id='1282', ...)]]

In this first filtering attempt we used a first order filter (input argument order=1). It can be seen, in the previous figure, that some unwanted information have been removed, unfortunately no filter has an ideal behavior, so despite we specify a high cutoff frequency of 40 Hz, some information above this threshold is maintained after filtering.

The good news are that components greater than 80 Hz are almost completely removed.

The filter performance can be improved by increasing the filter order, because the higher the filter order is, more quickly the transition between the pass and stop band will be. The transition band will be smaller because of a higher attenuation rate (-20 x order dB/decade).

However, the filter order must be chosen with precaution in order to avoid system instability. Magnitude Bode plots are very useful to check the filter response, as can be seen in the figure below, taking into consideration the following Mathematical formulation.

\begin{equation} G = -20\times\log_{10}\Bigg(\sqrt{1 + \bigg(\frac{f}{f_c}\bigg)^{2.n}}\Bigg) \end{equation}
$G$ - Gain factor (negative values reveal an attenuation) $n$ - Filter order (integer)
$f_{c}$ - Cutoff frequency of the filter (40 Hz, for the current implementation) $f$ - Independent variable (input frequency to be filtered)
In [11]: 
bsnb.plot_low_pass_filter_response()

7 - Repetition of the filtering stage but using a higher filter order

In [12]: 
# Digital low-pass filtering with a cutoff frequency f of 40 Hz
filter_signal_2 = bsnb.lowpass(signal, f=40, order=3, fs=sr)

# Power spectrum
freq_axis_3, power_spect_3 = bsnb.plotfft(filter_signal_2, sr)
In [13]: 
bsnb.plot_before_after_filter(signal, sr, band_begin=0.5, band_end=40, order=3, x_lim=[0, 100], y_lim=[0, 5e6],
                              orientation="same", show_plot=True)
Out[13]: 
[[Figure(id='2105', ...)]]

When the noise level is low, it may be difficult to observe its influence in time domain. In order to the digital filtering stage produce a visual effect in time domain, the noise level needs to be high.

So we will add some artificial noise do the signal and see the great impact of digital filtering.

E1 - Addition of artificial noise

In [14]: 
# Noise samples and translation of the baseline 
baseline = average(signal)
baseline_shift = 0.50 * baseline
noisy_signal = signal + normal(0, 1000, len(signal)) + baseline_shift

E2 - Noisy signal representation

In [15]: 
# Plotting of power spectrum    
bsnb.plot(linspace(0, len(noisy_signal) - 1, len(noisy_signal)), noisy_signal, x_axis_label='Sample Number', 
          y_axis_label='Raw Data', title="Noisy Signal", y_range=(-1e4, 6e4))

E3 - Digital Filtering Stage

In [16]: 
# Digital low-pass filtering with a cutoff frequency f of 40 Hz
noisy_signal_filter = bsnb.lowpass(noisy_signal, f=40, order=3, fs=sr)

E4 - Comparison of noisy and filtered signal in time domain

In [17]: 
bsnb.plot([linspace(0, len(noisy_signal) - 1, len(noisy_signal))]*2, [noisy_signal, noisy_signal_filter],
          grid_plot=True, grid_lines=1, grid_columns=2, x_axis_label='Sample Number', y_axis_label='Raw Data', 
          title=["Noisy Signal", "Filtered Signal"])

As described previously, none filter has an ideal behavior, but, in spite of not ideal the behavior of real filters is predictable.

For example, for the designed 3rd order Butterworth filter with cutoff frequency of 40 Hz, it is expected that after one decade (40 Hz x 10 = 400 Hz) the relative amplitude will be attenuated by -60 dB (assuming a value of 0.1% of the non-filtered signal).

In [18]: 
# Power spectrum (Noisy signal)
freq_axis_noisy, power_spect_noisy = bsnb.plotfft(noisy_signal, sr)

# Power spectrum (Filtered signal)
freq_axis_filter, power_spect_filter = bsnb.plotfft(noisy_signal_filter, sr)

# Relative amplitude 1 decade after the 40 Hz cutoff frequency --> 400 Hz
# Taking into consideration that the search is sequential, so, only the 
# first sample that meets the criterium is relevant.
index_decade = where(freq_axis_noisy >= 400)[0][0]
power_decade_noisy = power_spect_noisy[index_decade]
power_decade_filter = power_spect_filter[index_decade]
In [19]: 
from sty import fg, rs
print(fg(98,195,238) + "\033[1mRelative Amplitude/Power at 400 Hz frequency component [Noisy Signal]: \033[0m" + fg.rs + str(round(power_decade_noisy, 4)))
print(fg(148,193,30) + "\033[1mRelative Amplitude/Power at 400 Hz frequency component [Filtered Signal]: \033[0m" + fg.rs + 
      str(round(power_decade_filter, 4)))
print(fg(232,77,14) + "\033[1m\nRatio between filtered and noisy 400 Hz component [Attenuation]: \033[0m" + fg.rs + 
      str(round(power_decade_filter/power_decade_noisy, 4)) + " ~ " + str(round(power_decade_filter/power_decade_noisy, 3)) + " = 0.1 %")

Relative Amplitude/Power at 400 Hz frequency component [Noisy Signal]: 90071.9949
Relative Amplitude/Power at 400 Hz frequency component [Filtered Signal]: 80.2622

Ratio between filtered and noisy 400 Hz component [Attenuation]: 0.0009 ~ 0.001 = 0.1 %

Taking into consideration the previous demonstration, we can understand that the designed filter presents the desired behavior !

Unfortunately noise is everywhere and even physiological data can be classified as noise if we are studying signals with another nature, as a practical example, electromyographic (EMG) data will be noise when doing an electrocardiographic (ECG) acquisition.

But, as demonstrated before, we can face this obstacle with efficient solutions, using analogical (pre-acquisition) or digital (post-acquisition) filters. With this brief tutorial, it can be understood the basic functioning principle of digital filters and how the user should proceed to "design" a widely used Butterworth filtering system.

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

☌ Project Presentation
☌ GitHub Repository
☌ How to install biosignalsnotebooks Python package ?
☌ Signal Library 		☌ Notebook Categories
☌ Notebooks by Difficulty
☌ Notebooks by Signal Type
☌ Notebooks by Tag
In [20]: 
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()

.................... CSS Style Applied to Jupyter Notebook .........................
Out[20]: